• Tempo of Gene Expression
  • Build GCN (stepwise)
  • Build GCN (oneshot)
  • Rhy genes across tissues
  • Contact

On this page

  • ONE-SHOT (make_modules)
    • Make and plot modules
    • Modify network plot
  • Annotate the network
    • Identify rhythmic modules

Compare time-course gene expression data across tissues


Last updated on 2025-09-05.


Show code
library(dplyr)
library(dbplyr)
library(ggplot2)
for (i in list.files(here::here("R"), full.names = TRUE)) {
  source(i)
}

# SAMPLE NAME
## specify sample name
sample.names <- c(
  # dmel
  # "dmel-head",
  # "dmel-body",
  # "dmel-heart"
  # # mmus
  "mmus-brain_stem", 
  "mmus-cerebellum", 
  "mmus-hypothalamus", 
  "mmus-heart"
  # # panu
  # "panu-brain_stem",
  # "panu-cerebellum",
  # "panu-hypothalamus",
  # "panu-heart",
  # "panu-scn"
)
# sample.cycles <- c("LD", "DD")

## SPECIFY THE DATASET TO BUILD GCN WITH
which.sample <- sample.names[1]

writeLines(
  glue::glue("Reference tissue: {which.sample}")
)
Reference tissue: mmus-brain_stem
Structure of reference tissue data:

ONE-SHOT (make_modules)

Make and plot modules

Show code
# mods <- purrr::map(
#   sample.names,
#   ~ timecourseRnaseq::make_modules(
#       data = tidydata.db[[.x]],
#       log2 = TRUE,
#       id_column = "gene_name",
#       min_expression = NULL,  # automatically estimated
#       min_timepoints = NULL,  # automatically estimated
#       method = "wgcna",
#       qc = TRUE,
#       sim_method = "kendall",
#       soft_power = 15,      # automatically estimated
#       min_module_size = 50,
#       max_modules = 16,
#       merge_cutoff_similarity = 0.9,
#       plot_network_min_edge = 0.5,
#       plot_network = FALSE,
#       tidy_modules = TRUE
#     )
# ) |>
#   setNames(
#     sample.names
#   )

mods <- timecourseRnaseq::make_modules(
  data = tidydata.db[[which.sample]],
  log2 = TRUE,
  id_column = "gene_name",
  min_expression = NULL,  # automatically estimated
  min_timepoints = NULL,  # automatically estimated
  method = "wgcna",
  qc = TRUE,
  sim_method = "kendall",
  soft_power = NULL,      # automatically estimated
  min_module_size = 50,
  max_modules = 16,
  merge_cutoff_similarity = 0.9,
  plot_network = FALSE,
  # plot_network_min_edge = 0.5, # only used when `plot_network == TRUE`
  tidy_modules = TRUE
  )
---------------------------------------------------
1. Log2-transform and subset 
---------------------------------------------------
Applying log2-transformation...Done.
Estimated min_expression = 8 
Estimated min_timepoints = 6 
Subsetting data...Done.
[ NOTE ]: After subsetting, 8893 of 17406 rows remain.

 Flagging genes and samples with too many missing values...
  ..step 1
All okay!Visualizing the log-transformed data
---------------------------------------------------
2. Calculate similarity of expression 
---------------------------------------------------
Running gene-gene similarity...Done.


---------------------------------------------------
3. Create adjacency matrix 
---------------------------------------------------
Performing network topology analysis to pick
  soft-thresholding power...
   Power SFT.R.sq   slope truncated.R.sq mean.k. median.k. max.k.
1      1   0.7240  1.9500          0.654    4510    4900.0   6020
2      2   0.5090  0.5730          0.389    2970    3230.0   4790
3      3   0.0298  0.0544          0.363    2170    2280.0   4040
4      4   0.2950 -0.1720          0.710    1690    1680.0   3510
5      5   0.6060 -0.3070          0.785    1360    1280.0   3110
6      6   0.7420 -0.3970          0.817    1130     995.0   2800
7      7   0.8110 -0.4740          0.818     951     790.0   2540
8      8   0.8540 -0.5280          0.841     815     634.0   2330
9      9   0.8860 -0.5790          0.866     708     516.0   2140
10    10   0.8960 -0.6210          0.868     621     424.0   1980
11    12   0.9040 -0.6890          0.877     489     294.0   1730
12    15   0.9040 -0.7770          0.886     359     182.0   1440
13    18   0.8990 -0.8420          0.895     275     118.0   1220
14    21   0.8930 -0.8910          0.899     216      80.2   1050

Plotting resutls from the network topology analysis...

Done.
[ NOTE, FIGURE ]: Red horizontal line indicates a signed R^2 of 0.9

Setting soft-thresholding power to: 8.
Power-transforming the gene-gene similarity matrix...Done.

---------------------------------------------------
4. Convert into topological overlap matrix (dissTOM) 
---------------------------------------------------
Creating dissTOM...Done.
Performing hierarchical clustering on dissTOM...Done.

---------------------------------------------------
5. Identify modules (clusters) 
---------------------------------------------------
Merging modules that have a correlation ≥ 0.9 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
 antiquewhite4        bisque4           blue         coral1         coral2 
           215            669           1938             60            287 
          cyan    darkmagenta darkolivegreen    darkorange2        darkred 
           337            117            435             96            157 
 darkseagreen4  darkslateblue    floralwhite         grey60      honeydew1 
            61             87             97            172            388 
         ivory lavenderblush3      lightcyan     lightcyan1     lightgreen 
            98            155            174             98            165 
       magenta   mediumorchid  mediumpurple3   navajowhite2         orange 
           202             55            105             76            154 
    orangered4  paleturquoise palevioletred3           pink          plum1 
           108            123            301            623            110 
        purple            red         salmon        skyblue       skyblue2 
           199            234            186            127             54 
     steelblue          white        yellow4    yellowgreen 
           123            146             50            111 
Merging modules that have a correlation ≥ 0.85 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.

Module (cluster) size:
mergedColors
       bisque4         coral1       darkgrey    darkmagenta    darkorange2 
          1397            440           3195            117            194 
       darkred  darkseagreen4         grey60          ivory        magenta 
           157             61            172            672            325 
  mediumorchid  mediumpurple3   navajowhite2         orange palevioletred3 
           163            105             76            154            736 
         plum1    saddlebrown        skyblue       skyblue2          white 
           110            287            282             54            146 
       yellow4 
            50 
Merging modules that have a correlation ≥ 0.8 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Module (cluster) size:
mergedColors
      bisque4          blue        coral1   darkmagenta   darkorange2 
         2462          3195           440           796           194 
darkseagreen4         ivory       magenta  mediumorchid mediumpurple3 
           61           748           325           163           105 
       orange      skyblue2         white       yellow4 
          154            54           146            50 

Cutoff used: 0.8 
Number of modules identified: 14 

Calculating module-module similarity based
  on module-eigengene-expression...Done.
Tidying module names...Done.
Plotting adjacency matrix for module-module similarity...

---------------------------------------------------
6. Tidy modules (clusters) 
---------------------------------------------------

Modify network plot

Internal function; use ::: to call

Show code
# trash <- purrr::map(
#   which.sample,
#   ~ timecourseRnaseq:::plot_adj_as_network(
#       matrix = mods[[.x]][["adj_matrix_ME"]][["ME"]],
#       # layout = igraph::layout.sugiyama,
#       layout = igraph::layout_in_circle, # changed 
#       min_edge = 0.6,
#       node_label_size = 1.2,
#       node_size = 20,
#       edge_size = 3,
#       node_frame_col = "grey20",
#       node_fill_col = "grey80",
#       vertex.frame.width = 3
#     )
# )

timecourseRnaseq:::plot_adj_as_network(
  matrix = mods[["adj_matrix_ME"]][["ME"]],
  # layout = igraph::layout.sugiyama,
  layout = igraph::layout_in_circle, # changed 
  min_edge = 0.6,
  node_label_size = 1.2,
  node_size = 30,
  edge_size = 3,
  node_frame_col = "grey20",
  node_fill_col = "grey80",
  vertex.frame.width = 3
)
Visualizing a simplified representation of the circadian GCN

Annotate the network

Obtain a list of genes in each module

Show code
# module_genes <- purrr::map(
#   sample.names,
#   ~ mods[[.x]][["module_genes"]]
# ) |> 
#   setNames(sample.names)

module_genes <- mods[["module_genes"]]

Identify rhythmic modules

Show code
db_rhy <- purrr::map(
  sample.names,
  ~ load_rhy_genes(
      sample = .x
    )
) |> 
  setNames(sample.names)

# l_module_genes <- purrr::map(
#   sample.names,
#   ~ module_genes[[.x]] |> 
#     arrange(module_identity) |> 
#     group_split(module_identity) |> 
#     purrr::map(
#       ~ .x |> pull(gene_name)
#     ) |> 
#     setNames(unique(module_genes[[.x]]$module_identity))
# ) |> 
#   setNames(sample.names)

l_module_genes <- module_genes |> 
  arrange(module_identity) |> 
  group_split(module_identity) |> 
  purrr::map(
    ~ .x |> pull(gene_name)
  ) |> 
  setNames(unique(module_genes$module_identity))

###-###-###-###-###-###-###-###-
# Set your p-value of choice
# col_pval = "BH.Q"
col_pval = "default.pvalue"
# col_pval = "raw.pvalue"
###-###-###-###-###-###-###-###-

l_rhy_genes <- purrr::map(
  sample.names,
  ~ db_rhy[[.x]] |> 
    purrr::map(
      ~ .x |> 
        filter(
          if_all(
            all_of(col_pval),
            ~ .x < 0.05
          )
        ) |> 
        filter(
          ID %in% unlist(l_module_genes)
        ) |> 
        pull(1) |> 
        unique()
    ) |> 
    purrr::compact()
) |> 
  setNames(sample.names)

Modules vs. Rhythmic genes

Show code
# LIST ONE - WGCNA modules
list1 <- l_module_genes
sapply(list1, length) |> print()
  C1   C2   C3   C4   C5   C6   C7   C8   C9  C10  C11  C12  C13  C14 
 154  163 3195  105  748  146   54  796 2462  440  194   61   50  325 
Show code
trash <- purrr::map(
  sample.names,
  function(x) {
    cat("Tissue:", x, "\n")
    
    ## LIST TWO - rhythmic genes
    list2 <- l_rhy_genes[[x]]
    sapply(list2, length) |> print()
    
    ## CHECK FOR OVERLAP
    # define size of genome
    size = length(unique(c(unlist(list1), unlist(list2))))
    # make a GOM object
    gom.1v2 <- GeneOverlap::newGOM(
      list2, 
      list1,
      genome.size = size
    )
    GeneOverlap::drawHeatmap(
      gom.1v2,
      adj.p = TRUE,
      cutoff = 0.05,
      what = "odds.ratio",
      # what="Jaccard",
      log.scale = T,
      note.col = "black",
      grid.col = "Oranges"
    )
    
    gom.1v2
  }
)
Tissue: mmus-brain_stem 
      ARS    empJTK GeneCycle       JTK    meta2d      RAIN 
     1010       361       368       198       783       687 

Tissue: mmus-cerebellum 
      ARS    empJTK GeneCycle       JTK    meta2d      RAIN 
     2663       701       641       241      1952       593 

Tissue: mmus-hypothalamus 
      ARS    empJTK GeneCycle       JTK    meta2d      RAIN 
     2042       559       481       223      1294       609 

Tissue: mmus-heart 
      ARS    empJTK GeneCycle       JTK    meta2d      RAIN 
     2550       987       785       532      2251      1316